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Abstract 

Phonon plays essential roles in dynamical behaviors and thermal properties, which are central topics in fundamental 
issues of materials science. The importance of first principles phonon calculations cannot be overly emphasized. Phonopy 
is an open source code for such calculations launched by the present authors, which has been world-widely used. Here we 
demonstrate phonon properties with fundamental equations and show examples how the phonon calculations are applied 
in materials science. 


1. Introduction 

Application of first principles calculations in condensed 
matter physics and materials science has greatly expanded 
when phonon calculations became routine in the last decade 
Thanks to the progress of high performance computers and 
development of accurate and efficient density functional 
theory (DFT) codes, a large set of first principles calcu¬ 
lations are now practical with the accuracy comparable 
to experiments using ordinary PC clusters. In addition 
to electronic structure information, a DFT calculation for 
solids provides energy and stress of the system as well as 
the force on each atom. Equilibrium crystal structures can 
be obtained by minimizing residual forces and optimizing 
stress tensors. When an atom in a crystal is displaced 
from its equilibrium position, the forces on all atoms in 
the crystal raise. Analysis of the forces associated with a 
systematic set of displacements provides a series of phonon 
frequencies. First principles phonon calculations with a fi¬ 
nite displacement method can be made in this way. The 
present authors have launched a robust and easy-to-use 
open-source code for first principles phonon calculations, 

phonopy [lliaialllligEllzllHlEllinillllliailsllIllllSlIle]. 

The number of users is rapidly growing world-wide, since 
the information of phonon is very useful for accounting 
variety of properties and behavior of crystalline materials, 
such as thermal properties, mechanical properties, phase 
transition, and superconductivity. In this article, we show 
examples of applications of the first principles phonon cal¬ 
culations. 

In Sections [2]|^ we take FCC-Al as examples of ap¬ 
plications of first principles phonon calculations. For the 
electronic structure calculations, we employed the plane- 
wave basis projector augmented wave method m in the 
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framework of DFT within the generalized gradient approx¬ 
imation in the Perdew-Burke-Ernzerhof form |18j as im¬ 
plemented in the VASP code [HI EQI [2T] . A plane-wave 
energy cutoff of 300 eV and an energy convergence criteria 
of 10“® eV were used. A 30 x 30 x 30 fc-point sampling 
mesh was used for the unit cell and the equivalent density 
mesh was used for the supercells together with a 0.2 eV 
smearing width of the Methfessel-Paxton scheme [22]. Eor 
the phonon calculations, supercell and finite displacement 
approaches were used with 3x3x3 supercell of the conven¬ 
tional unit cell (108 atoms) and the atomic displacement 
distance of 0.01 A. 

2. Harmonic approximation 

In crystals, it is presumed that atoms move around 
their equilibrium positions r{lK) with displacements u((k), 
where I and n are the labels of unit cells and atoms in 
each unit cell, respectively. Crystal potential energy $ is 
presumed to be an analytic function of the displacements 
of the atoms, and $ is expanded as 

$ = ^>0 + X X ^c{ln)ua{lK) 

Ik oc 

+ ^ X I'k')uc{Ik)up{ 1'n') 

W kk' aj3 



ll'l" kk'k" Oi^-y 

I' k' , I” K!')Ua{lK)ui}{l' k”) -(-••• ( 1 ) 

where a, j3, ■ ■ ■ are the Cartesian indices. The coeffi¬ 
cients of the series expansion, Z'/c'), 

and, ^ai 3 -y{lK,l' k',1”k”), are the zero-th, first, second, and 
third order force constants, respectively. With small dis¬ 
placements at constant volume, the problem of atomic 
vibrations is solved with the second-order terms as the 
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harmonic approximation, and the higher order terms are 
treated by the perturbation theory. 

With a force Fq(/k) = — an element of second- 

order force constants k') is obtained by 




dua{ln)du^{l' K,’) 


dFpjl’K’) 

dUailK) 


Crystal symmetry is utilized to improve the numerial ac¬ 
curacy of the force constants and to reduce the computa¬ 
tional cost. The more details on the calculation of force 
constants are found in Ref. HE]. 

As found in text books [231 [211 ESI ES], dynamical prop¬ 
erty of atoms in the harmonic approximation is obtained 
by solving eigenvalue problem of dynamical matrix D(q), 


D(q)eqj = e, 




or 
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where is the mass of the atom k, q is the wave vector, 
and j is the band index. lOqj and eqj give the phonon fre¬ 
quency and polarization vector of the phonon mode labeled 
by a set {q, j}, respectively. Since D(q) is an Hermitian 
matrix, its eigenvalues, w^-, are real. Usually D(q) is ar¬ 
ranged to be a Sria x 3na matrix |26j . where 3 comes from 
the freedom of the Cartesian indices for crystal and Ua is 
the nubmer of atoms in a unit cell. Then Gq^ becomes a 
complex column vector with 3na elements, and usually Gq^ 
is normalized to be 1, i.e., 1®^!^ = 1- ®qi contains 

information of collective motion of atoms. This may be 
understood as a set of atomic displacement vectors. 


{u(n),... ,u(/k) } = 

/ ^ pi p*qr(ii) _ 


A 


ria iqr(iK) 


(5) 


where A is the complex constant undetermined by Eq. ([^ , 
and G^ A' = p^^ p^^l 

diiu Cqj v'^qj ’ '"qj ’ ' ’ 

As a typical example, the phonon band structure and 
phonon density of states (DOS) of A1 are shown in Fig. 
The phonon DOS is defined as 


5(w) 



^qj) I 


( 6 ) 


where N is the number of unit cells in crystal. Divided by 
N^ g{uj) is normalized so that the integral over frequency 
becomes Sria. The phonon band structure can be directly 
comparable with experimental data by neutron or X-ray 
inelastic scattering. They often show reasonable agree¬ 
ments [HEZIEH]- Frequency data by Raman and infrared 
(IR) spectroscopy can also be well reproduced (5] ES] ■ Ir¬ 
reducible representations of phonon modes, which can be 



Figure 1: (color online) Phonon band structure and DOS of Al. 


used to assign Raman or IR active modes, are calculated 
from polarization vectors li ES]. Atom specific phonon 
DOS projected along a unit direction vector n is defined 
as 

5K(w,n) = ^^(5(a;-a;qj-)|h.G^-|2. (7) 

qj 

This gK,{u},n) can be directly compared with that mea¬ 
sured by means of nuclear-resonant inelastic scattering us¬ 
ing synchrotron radiation. In Ref. m, phonon calcula¬ 
tions of Llg-type FePt projected along the c-axis and basal 
plane are well comparable to experimental ®^Fe nuclear- 
resonant inelastic scattering spectra measured at 10 K in 
the parallel and perpendicular geometries, respectively. 



Figure 2: (color online) Thermal properties of Al. Entropy, Cy, and 
Helmholtz free energy were calculated with harmonic approximation 
(Sec.j^. QHA was employed to obtain Cp (Sec.j^. Physical units 
are shown with labels of the physical properties, and the value of the 
vertical axis is shared by them. Dotted curve depicts the experiment 
of Cp [31]. 


Once phonon frequencies over Brillouin zone are known, 
from the canonical distribution in statistical mechanics for 
phonons under the harmonic approximation, the energy E 
of phonon system is given as 


E — hujqj 

qj 


1 

2 


1 

exp(^a;q,/kBT) - 1 


( 8 ) 
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where T, k^, and Ti are the temperature, the Boltzmann 
constant, and the reduced Planck constant, respectively. 
Using the thermodynamic relations, a number of ther¬ 
mal properties, such as constant volume heat capacity Cy, 
Helmholtz free energy F, and entropy S, can be computed 
as functions of temperature [55]: 

exp(;ta;qj/fcBT) 

(9) 

F = ^ ^ hLOqj + knT ^ ln[l - exp(-?ia;qj7fcBr)], 
qi qj 

( 10 ) 

and 

^ [;ia;qj/2fcBr] 

q? 

— /cb In [2 sinh(?ta;qj /2ky,T)\ . (11) 

qi 

The calculated F, Cy, and S for A1 are shown in Fig.j^ 


4. Quasi-harmonic approximation 


By changing volume, phonon properties vary since the 
crystal potential is an anharmonic function of volume. 
In this article, the term “quasi-harmonic approximation 
(QHA) “ means this volume dependence of phonon proper¬ 
ties, but the harmonic approximation is simply applied at 
each volume. Figure]^ shows calculated phonon frequen¬ 
cies of A1 at X and L points with respect to ten different 
unit-cell volumes. Typically phonon frequency decreases 
by increasing volume, and the slope of each phonon mode 
is nearly constant in the wide volume range. The nor¬ 
malized slope is called mode-Griineisen parameter that is 
defined as 


1^{V) 


V dw^,{V) 

dV 


(15) 


Once dynamical matrix is known, 7 qj is easily calculated 
from the volume derivative of the dynamical matrix [55] . 
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3. Mean square atomic displacements 


With the phase factor convention of the dynamical ma¬ 
trix used in Eq. Q, an atomic displacement operator is 
written as. 


Uailn) 
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where and d are the creation and annihilation oper¬ 
ators, respectively. Distributions of atoms around their 
equilibrium positions are then obtained as the expectation 


values of Eq. (12|. The mean square atomic displacement 


projected along h is obtained as 


(|Ufi(/c)7) = 


1 + 277^0 


2Nm^ ^ 

q? 


n • e” 




(13) 


Eq. (131 is lattice-point (1) independent since the phase 


factor disappears. Anisotropic atomic displacement pa¬ 
rameters (ADPs) to estimate the atom positions during 
thermal motion can also be computed and compared with 
experimental neutron diffraction data. Thermal ellipsoids 
may be discussed using mean square displacement matrix 
B(/c) defined by 


B(«;) = 


2Nmn 


1 + 2nqj- 


UJa 


■'qj ■ 


(14) 


The shape and orientation of an ellipsoid is obtained solv¬ 
ing eigenvalue problem of this matrix. The method has 
been applied to show the ORTEP (Oak Ridge Thermal 
Ellipsoid Plot)-style drawing of ADPs [14]. Ref. [7] shows 
an example for a ternary carbide Ti 3 SiC 2 having a layered 
structure known as MAX phases, in which we can see good 
agreement between calculated and experimental ADPs. 


The quantity can be related to macroscopic Griineisen pa¬ 
rameter 7 using mode contributions to the heat capacity 
Cqj found in Eq. 0,7 = X)qj Im^qjlCy [511 [55]. 

Silicon is known as a famous exception to have large 
negative mode-Griineisen parameters. Mode-Griineisen pa¬ 
rameter is a measure of anharmonicity of phonon modes 
and is related to third-order force constants directly [55] . 
Therefore crystals that possess large anharmonic terms 
beyond third-order terms in Eq. Q can show non-linear 
change of phonon frequency with respect to volume. This 
is often observed in crystals that exhibit second- or higher- 
order structural phase transitions [5]. 

The phonon frequency influences the phonon energy 
through Eq. 0. The thermal properties are thereby af¬ 
fected. Using thermodynamics definition, thermodynamic 
variables at constant volume is transformed to those at 
constant pressure that is often more easily measurable in 
experiments. Gibbs free energy G{T,p) at given temper¬ 
ature T and pressure p is obtained from Helmholtz free 
energy F(T; V) through the transformation, 

G(F,p)=mm[F(F;l/)+pU], (17) 

where the right hand side of this equation means finding a 
minimum value in the square bracket by changing volume 
V. We may approximate F(F; V) by the sum of electronic 
internal energy Uei{V) and phonon Helmholtz free energy 
Fph(T; U), i.e., F{T-, U) ~ [/ei(U) + Fph(T; U). UeiiV) is 
obtained as total energy of electronic structure from the 
first principles calculation, and the first principles phonon 
calculation at T and V gives Fph(r; V). The calculated 
Ue\{V) -I- Fph(T; V) are depicted by the filled circle sym¬ 
bols in Fig. Wp, where the ten volume points chosen are 
the same as those in Fig. The nine curves are the 
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Figure 3: (color online) (a) Phonon frequencies of A1 at X and L 
points with respect to unit cell volume are shown by filled and open 
circles, respectively. The solid and dotted lines are guides to the eye. 

(b) Uel + -Fph of A1 with respect to volume at temperatures from 
0 to 800 K with 100 K step are depicted by filled circles and the 
values are fit by the solid curves. Cross symbols show the energy 
bottoms of the respective curves and simultaneously the equilibrium 
volumes. Lines connecting the cross symbols are guides to the eye. 

(c) Volumetric thermal expansion coefficient of Al. Calculation is 
shown with solid curve and experiments are depicted by filled circle 
symbols m and dotted curve m- 


fits to equation of states (EOS) at temperatures from 0 
to 800 K with 100 K step. Here the Vinet EOS [34] was 
used to fit the points to the curves though any other func¬ 
tions can be used for the fitting. The minimum values 
at the temperatures are depicted by the cross symbols in 
Fig.[3)o and are the Gibbs free energies at the tempera¬ 
tures and the respective equilibrium volumes are simulta¬ 
neously given. Volumetric thermal expansion coefficient, 
/^(F) = v(V) ^ obtained from the calculated equi¬ 

librium volumes V(T) at dense temperature points. I3(T) 
for Al is shown in Fig.j^, where we can see that the calcu¬ 
lation shows reasonable agreements with the experiments. 
In thermodynamics, heat capacity at constant pressure Cp 
is given by 


Cp{T,p) = -T 


d^G{T,p) 


dT2 


= Cv{T,ViT,p))+T 


dV{T,p) dS{T- V) 


dT 


dV 


V=V(T,p) 


■ ( 18 ) 


In Eq. (18), the second term of the second equation is un¬ 


derstood as the contribution to heat capacity from thermal 
expansion. Cp for Al is shown in Eig. The agreement 
of the calculation with the experiment is excellent. At 
high temperatures, the difference between Cp and Cy is 
not negligible in Al. Therefore it is essential to consider 
thermal expansion for heat capacity. 

QHA is known as a reasonable approximation in a wide 
temperature range below melting point except for temper¬ 
atures very close to melting point where higher-order terms 
in Eq. 0 become evident [55] . 


5. Stability condition and imaginary mode 

At equilibrium, = 0, a crystal is dynamically 

(mechanically) stable if its potential energy always increases 
against any combinations of atomic displacements. In the 
harmonic approximation, this is equivalent to the condi¬ 
tion that all phonons have real and positive frequencies [55] . 
However under virtual thermodynamic conditions, imag¬ 
inary frequency or negative eigenvalue can appear in the 
solution of Eq. 0. This indicates dynamical instability of 
the system, which means that the corrective atomic dis¬ 
placements of Eq. ([^ reduce the potential energy in the 
vicinity of the equilibrium atomic positions. 


(a) a-phase (FCC) (b) (J-phase (BCC) 



(c) m-phase 



Figure 4: (color online) Phonon band structures of (a) a-Ti (HCP), 
(b) /?-Ti (BCC), and (c) a;-Ti. The figure (d) shows the hexagonal 
crystal structure of ta-Ti. 


Imaginary mode provides useful information to study 
displacive phase transition. A typical example is shown 
in Figs. to [12] ■ Imaginary modes can be found 

only for /?-Ti, that has BCC structure, at both P and N 
points. This indicates that /3-Ti is unstable at low tem¬ 
perature. Such imaginary modes cannot be seen for either 
w-Ti whose crystal structure is shown in Fig. [^ or a- 
Ti (HCP). Experimentally /3-Ti is known to occur above 
1155K. At such high temperatures, large atomic displace¬ 
ments can stabilize the BCC structure. In such a case, 
the perturbation approach is invalid. Phonons with large 
atomic displacements may be treated by self-consistent 
phonon method [551 EH] or by a combination of molecu¬ 
lar dynamics and lattice dynamics calculation [37l 138] [39] . 
which is not discussed in this article. 

A given structure having imaginary phonon modes can 
be led to alternative structures through continuous atomic 
displacements and lattice deformations. The present au¬ 
thors systematically investigated the evolution of crystal 
structures from the simple cubic (SC) structure [6]. The 
inset of Fig. [^ shows the phonon band structure of SC-Cu 
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{Pm3m). Imaginary modes can be found at M(l/2,1/2,0) 
and X(l/2,0,0) points. Then the SC structure is de¬ 
formed along these directions. In order to accommodate 
the deformation in the calculation model with the peri¬ 
odic boundary condition, the unit cells are expanded by 
2 X 2 X 1 for the M point and 2 x 1 x 1 for the X point. The 
M point deformation breaks the crystal symmetry of SC 
(PmSm) to PA/nmm. The doubly degenerated instabil¬ 
ity at the X point leads to Pmma and Cmcm as highest 
possible symmetries. The deformed crystal structures are 
relaxed by first principles calculations imposing the corre¬ 
sponding space-group operations. After these procedures, 
body-centered tetragonal (BCT), simple hexagonal (SH), 
and FCC are respectively formed. The whole procedure 
finishes when all crystal structures at the end-points are 
found to be dynamically stable. Finally a treelike struc¬ 
ture of crystal structure relationships was drawn as shown 
in Fig. where thick lines indicate phase transition path¬ 
ways (PTPs). The space-group type written near a line 
is a common subgroup of initial and final structures. The 
presence of the line indicates that the energy decreases 
monotonically with the phase transition. In other words, 
the transition can take place spontaneously without any 
energy barrier. The line is terminated when the final struc¬ 
ture is dynamically stable. 



Figure 5: (color online) Line diagram of structural transition path¬ 
ways in Cu. The inset shows phonon band structure of simple cubic 
(SC) Cu. Open and filled symbols represent dynamically unstable 
and stable crystal structures, respectively. Lines connecting these 
symbols are the phase transition pathways for which space-group 
types are shown near the lines. 

In the line diagram, uj is located at the junction of two 
pathways, i.e., w BCC HCP and lo FCC. The in¬ 
stability of w at the F point leads to BCC, which is still dy¬ 
namically unstable and eventually leads to HCP. Another 
instability at the M point leads to FCC. The other insta¬ 
bility at the K point, which is doubly degenerate, leads to 
FCC. On the path from ui to BCC, the crystal symmetry 
of UJ having the space-group type of P6/mmm is once low¬ 
ered to P3ml and then becomes Im3m (BCC) after the 


geometry optimization. Both uj- and BCC-Cu are dynam¬ 
ically unstable, which can be formed only under crystal 
symmetry constraints. FCC-Cu is, of course, dynamically 
stable. Several PTPs between BCC and FCC have been 
proposed in literature. However, they are mostly based 
only upon investigation of continuous lattice deformation. 
For example in the classical Bain path, formation of BCT 
in between BCC and FCC is considered. Formation of SC 
is taken into account in the “trigonal Bain path.” Normal 
modes of phonon, which should be most adequate to de¬ 
scribe the collective atomic displacements, have not been 
considered. The presence of uj as the lowest energy barrier 
in the BCC-FCC pathway had not been reported before 
Ref. [5]. The situation is the same for the BCC-HCP tran¬ 
sition, known as the Burgers path. The Burgers path was 
thought to be quite complicated from the viewpoint of the 
lattice continuity. On the basis of the present study, it 
can be easily pointed out that the BCC-HCP transition 
pathway is along the space-group type of Cmcm. 

Evolution diagram was constructed in the same way for 
Nai?Ti 04 {R: rare-earth metal) with Ruddlesden-Popper 
type structure [9] . Inversion symmetry breaking by oxygen 
octahedral rotations was unambiguously demonstrated. The 
mechanism is expected to lead to many more families of 
acentric oxides. 

6. Interactions among phonons and lattice thermal 

conductivity 

Using the harmonic phonon coordinates, anharmonic 
terms in Eq. 0 are transformed to a picture of phonon- 
phonon interactions Lattice thermal conductivity 

can be accurately calculated by solving linearized Boltz¬ 
mann transport equation with the phonon-phonon interac¬ 
tion strength obtained using first principles calculation [3 
SOI SI]. Although the computational cost for such calcula¬ 
tions is many orders of magnitudes higher than the ordi¬ 
nary DFT calculations of primitive cells, such calculations 
have already been applied for many simple compounds and 
computed lattice thermal conductivities show good agree¬ 
ments with experimental data 0 m]. Calculations with 
special focus on searching thermoelectric materials have 
also been made HnHH Ell- 
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